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Abstract 

We  present  rigorous  a  posteriori  error  bounds  for  the  Empirical  Interpolation  Method  (EIM).  The  essential 
ingredients  are  (i)  analytical  upper  bounds  for  the  parametric  derivatives  of  the  function  to  be  approximated,  (ii) 
the  EIM  “Lebesgue  constant,”  and  (in)  information  concerning  the  EIM  approximation  error  at  a  finite  set  of 
points  in  parameter  space.  The  bound  is  computed  “offline”  and  is  valid  over  the  entire  parameter  domain;  it  is 
thus  readily  employed  in  (say)  the  “online”  reduced  basis  context.  We  present  numerical  results  that  confirm  the 
validity  of  our  approach. 

Resume 

Un  estimateur  a  posteriori  d’erreur  pour  la  methode  d’interpolation  empirique.  On  introduit  des 
bornes  d’erreur  a  posteriori  rigoureuses  pour  la  methode  d’interpolation  empirique,  EIM  en  abrege  (pour  Empirical 
Interpolation  Method).  Les  ingredients  essentiels  sont  (i)  des  bornes  analytiques  des  derivees  par  rapport  au 
parametre  de  la  fonction  a  interpoler,  (ii)  une  “constante  de  Lebesgue”  de  EIM,  et  (in)  de  l’information  sur  l’erreur 
d’approximation  commise  par  EIM  en  un  nombre  fini  de  points  dans  l’espace  des  parametres.  La  borne,  une  fois 
pre-calculee  “hors-ligne” ,  est  valable  sur  tout  l’espace  des  parametres ;  elle  peut  done  etre  utilisee  directement 
telle  quelle  dans  les  applications  (etape  “en  ligne”  des  calculs  dans  le  contexte  de  la  methode  des  bases  reduites). 
On  montre  des  resultats  numeriques  qui  confirment  la  validite  de  notre  approche. 
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Version  frangaise  abregee 

Soit  p)  =  Yl!k= 1  4>m{p)F{-\  Pm)  une  approximation  cle  T  £  C°°(T>,  fonction  qui  joue  le 

role  d’un  coefficient  parametre  “non-affine”  dans  la  methode  des  bases  reduites.  La  methode  d’interpola- 
tion  empirique  (EIM)  sert  a  construire  p)  j  elle  fournit  aussi  un  estimateur  de  l’erreur  d’interpolation 

ej v/(m)  —  —  J7m(s  M)IU°°(n))  qui  n’est  pas  rigoureux  mais  qui  est  souvent  suffisamment  precis. 

Dans  ce  travail,  nous  construisons  rigoureusement  une  borne  superieure  d’erreur  a  posteriori  pour  cm(m)- 

D’abord,  nous  rappelons  ce  qu’est  la  methode  EIM.  Ses  ingredients  essentiels  sont  (i)  la  construc¬ 
tion  d’un  espace  d’approximation  Wm  =  span{jF(-;^m)}^=1  avec  quelques  valeurs  pm,  1  <  m  <  M, 
selectionnees  par  un  algorithme  glouton,  pour  le  parametre  p  £  T>,  et  (ii)  la  selection  d’un  ensemble  de 
noeuds  d’interpolation  Tm  =  {ti  £  fi, . . . ,  4m  €  Q}  associe  a  Wm .  L’approximation  Tm ( ' ;  p)  £  Wm  est 
definie  comme  l’interpolant  de  p)  sur  l’ensemble  Tm- 

Ensuite,  nous  introduisons  notre  nouvelle  borne  d’erreur.  Pour  cela,  nous  developpons  T[x\  p)  en  une 
serie  de  Taylor  a  plusieurs  variables,  avec  un  ensemble  fini  <I>  de  points  dans  l’espace  des  parametres 
V  C  Rp.  Pour  tout  entier  positif  I,  nous  supposons  max^gD  max^^p  ||<F^(-;M)IU00(n)  <  oy  (<  oo), 
avec  Mf  l’ensemble  de  tous  les  multi-indices  positifs  (3  =  (fii, . . .  f3p)  de  dimension  P  et  de  longueur 
J^iLi  Pi  —  I  (pour  1  <  i  <  P,  Pi  est  un  entier  positif)  et  la  derivee  /3-ieme  de  T  par  rapport  a 
p.  Puis  nous  posons  =  max^g-p  minr6$  \p  —  r|,  nous  clefinissons  une  “constante  de  Lebesgue”  A m  = 
suPcceaX)m=i  l^mf(a;)l>  °il  Vm(.x)  €  Wm  sont  les  fonctions  caracteristiques  V^(tn)  =  6mn,  1  <  m,n  < 
M,  et  nous  pouvons  alors  prouver  notre  Proposition  2.1,  soit  :  ma <  Sm,p,  avec  une  borne 
Sm,p  definie  en  (2). 

Enfin,  nous  presentons  des  resultats  numeriques  avec  une  fonction  gaussienne  T{x\  p)  =  exp(—  ((xi  — 
Xi)2  +  (x2  —  X2)2) /% <a2)  sur  O  =  (0,  l)2.  Avec  un  seul  parametre  (scalaire)  a  =  p  £  T>\  =  [0.1,1]  et 
(xi,x2)  —  (0.5, 0.5)  fixe,  on  calcule  max^gs^  eM(p)  et  6m, p,  P  =  1, 2, 3, 4  pour  1  <  M  <  Mmax  (Fig.  1). 
Nous  observons  que  les  bornes  d’erreur  commencent  par  decroitre  puis  atteignent  un  “plateau”  en  M . 
On  calcule  aussi  A m  et  l’effectivite  moyenne  f}M,p  pour  p  =  4  en  tant  que  fonctions  de  M  :  la  constante 
de  Lebesgue  ne  croit  que  legerement  avec  M,  et  les  bornes  d’erreurs  sont  tres  precises  pour  de  petites 
valeurs  de  M.  Avec  deux  parametres  {x\,X2)  =  p  £  T>n  =  [0.4, 0.6]2  et  a  =  0.1  fixe,  les  resultats  sont 
similaires  au  cas  d’un  seul  parametre  (Fig.  2). 


1.  Introduction 

The  Empirical  Interpolation  Method  (EIM),  introduced  in  [1],  serves  to  construct  “affine”  approxima¬ 
tions  of  “non-affine”  parametrized  functions.  The  method  is  frequently  applied  in  reduced  basis  approxi¬ 
mation  of  parametrized  partial  differential  equations  with  non-affine  parameter  dependence  [4] ;  the  affine 
approximation  of  the  coefficient  functions  is  crucial  for  computational  efficiency.  In  previous  work  [1,4]  an 
estimator  for  the  interpolation  error  is  developed;  this  estimator  is  often  very  accurate,  however  it  is  not  a 
rigorous  upper  bound.  In  this  paper,  we  develop  a  rigorous  a  posteriori  upper  bound  for  the  interpolation 
error  and  we  present  numerical  results  that  confirm  the  validity  of  our  approach. 

To  begin,  we  summarize  the  EIM  [1,4].  We  are  given  a  function  Q  :  Q  x  V  — *  R.  such  that,  for  all  p  £  V, 
Q(-:  p)  £  here,  T>  C  is  the  parameter  domain,  Vt  C  R2  is  the  spatial  domain — a  point  in  which 

shall  be  denoted  by  x  =  {x\,X2) — and  L°°(Cl)  =  {v  |esssup„eQ  |w(x)|  <  00}.  We  introduce  a  finite  train 
sample  5trajn  C  V  which  shall  serve  as  our  V  surrogate,  and  a  triangulation  T/y(fl)  of  fi  with  J\f  vertices 
over  which  we  shall  in  practice  realize  Q(-;  p)  as  a  piecewise  linear  function. 

We  first  define  the  nested  EIM  approximation  spaces  Wfj,  1  <  M  <  Mmax.  We  first  choose  pi  €  V, 
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compute  gi  =  G(-',g  1),  and  define  Wf  =  spanjgi};  then,  for  2  <  M  <  Mmax,  we  determine  gM  = 
arg maxMeHtrain  in tzeWe_i  ||£(-;  b)  -  compute  gM  =  G{ ■;  Mm)j  and  define  =  span{gm}^f=1. 

We  next  introduce  the  nested  set  of  EIM  interpolation  nodes  =  {t\, . . .  ,tM},  1  <  M  <  Mmax. 
We  first  set  t\  =  argsup^.gp  | r/i (a;) |  and  qi  =  gi/gi{t\)\  then,  for  2  <  M  <  Mmax,  we  solve  the  linear 
system  =  9m(L),  1  <  i  <  M  -  1,  and  set  rM(x)  =  gM{x )  -  X^li1  w )>  = 

argsup^gQ  |rM(a;)|,  and  <7m  =  For  1  <  M  <  Mmax,  we  define  the  matrix  B M  G  RMxM 

such  that  B™  =  qj(ti),  1  <  i,j  <  M;  we  note  that  BM  is  lower  triangular  with  unity  diagonal  and  that 

Um}m=  l  is  a  basis  for  I1.4]- 

We  are  now  given  a  function  TL  :  O  x  V  — >  R.  such  that,  for  all  g  £  V,  TL{-\g)  €  L°°(Q).  We  define 
for  any  g  £  V  the  EIM  interpolant  Ttwer(-;g)  £  Wf{  as  the  interpolant  of  H(-\ g)  over  the  set  T^- 

Specifically  'Hws[ (saO  =  where  J2jLi  Bij  j(/i)  =  Hitgg),  1  <i<  M.  Note  that 

the  determination  of  the  coefficients  </>Mm(p)  requires  only  0{M2)  computational  cost. 

Finally,  we  define  a  “Lebesgue  constant”  [6]  A m  =  sup,,,  efi  Em=l  \Vm{x)l  where  Vm  €  WM  are  the 
characteristic  functions  of  Wfr  satisfying  V^f  ( tn )  =  Smn,  1  <  to,  n  <  M;  here,  Smn  is  the  Kronecker  delta 
symbol.  We  recall  that  (i)  the  set  of  all  characteristic  functions  is  a  basis  for  Wfr ,  and  (ii) 

the  Lebesgue  constant  A m  satisfies  A m  <  2M  —  1  [1,4].  In  applications,  the  actual  asymptotic  behavior 
of  Am  is  much  better,  as  we  shall  observe  subsequently. 


2.  A  Posteriori  Error  Estimation 

We  now  develop  the  new  and  rigorous  upper  bound  for  the  error  associated  with  the  empirical  inter¬ 
polation  of  a  function  T  :  x  V  — »  M.  We  shall  assume  that  T  is  parametrically  smooth;  for  simplicity 
here,  we  suppose  T  £  C°°(T>,  L°°(Gl)).  Our  bound  depends  on  the  parametric  derivatives  of  T  and  on 
the  EIM  interpolant  of  these  derivatives.  For  this  reason,  we  introduce  a  multi-index  of  dimension  P, 
(3  =  /?p),  where  the  /3j,  1  <  i  <  P,  are  non-negative  integers;  we  further  define  the  length 

\/3\  =  Eti  A,  and  denote  the  set  of  all  distinct  multi-indices  /3  of  dimension  P  of  length  /  by  Mf .  The 
cardinality  of  Mf  is  given  by  card(A if)  =  (P+/_1).  For  any  multi-index  /3,  we  define 

f)\0\  T 

Fi0)(x;g)  =  — g- - (1) 

^(i)  •  •  •  M(P) 

we  require  that  max^gp  max^g^p  ||pW(-;  /z)||L°°(n)  <  crp  (<  °°)  for  non-negative  integer  p. 

Given  any  g  £  V,  we  define  for  1  <  M  <  Mmax  the  interpolants  of  P(-;  g)  and  g)  as  Pm(-;  f-)  — 

Pwr(--,g)  and  (P^)m(-;  g)  =  pi^(-;/z),  respectively.  We  emphasize  that  both  interpolants 

M  WM 

and  (P^))m( Ll)  lie  in  the  same  space  W^j — we  do  not  introduce  a  separate  space,  ,  spanned 

by  solutions  of  1  <  M  <  Mmax.  It  is  thus  readily  demonstrated  that  0)  = 

A*),  which  we  thus  henceforth  denote  Pj^(-; /z). 1  Note  that  (•;  g)  £  W is  the  unique 
interpolant  satisfying  =  P^l(tTO;/z),  1  <  rn  <  M.  We  can  further  demonstrate  [3]  in  certain 

cases  that  if  Tm{-\ g)  tends  to  P(-;  g)  as  M  — >  00  then  P^(-;  g)  tends  to  g)  as  M  — >  00. 


1  Let  Zq  =  [<ji-..<2m]  and  tM  =  [ti...tM].  We  then  have  FM(-\A)  =  Zq(BM)  1F(tM\A)  and  A)  = 

Since  BM  and  the  basis  functions  <{, .  1  <  i  <  M,  are  independent  of  //.  it  follows  that 

(^m)(/3)( -;m)  =  {Zq{BM)-'T{iM-,p))W  =  ZqiBM)-1^)^^)  = 
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We  now  develop  the  interpolation  error  upper  bound.  To  begin,  we  introduce  a  set  of  points  d>  C  T>  of 
size  and  define  p$  =  max^gp  minre$  \p  —  r|;  here  |  •  |  is  the  usual  Euclidean  norm.  We  then  define 

8m, P  =  (1  +  A M)^j  pI  Pp/ 2  +  sup  (  PJ/ 2  maXP  H-F(/?)(-; T)  _  P<m^ T)l U“(0)  ]  •  (2) 

P-  re*  \j=o  J'  peM<  ) 

We  can  now  demonstrate 

Proposition  2.1  For  given  positive  integer  p,  max^-p  ||P(-;  p)  —  JFm{-\  p)\\ <  8m, p,  V/jeP,  1  < 


M  <  M„ 


PROOF.  We  present  the  proof  for  P  =  1  and  refer  the  reader  to  [3]  for  the  general  case  P  >  1.  For 
brevity,  we  first  define  (assuming  existence)  Ag(T,p)  =  Y^'jZ o  r)  —  as  the  first  p  terms  in  the 

Taylor  series  of  Q  around  r.  We  then  choose  r  as  r*(p)  =  argmin?6$  \p  —  f  |.  We  note  that 

11^0; p)  - ^M(-;^)IL-o(n)  <  m)IU(q)  +  \\apAt*,p)-  ^m(-;m)IILoo(0)  (3) 

for  all  p  G  T>.  We  recall  the  univariate  Taylor  series  expansion  with  remainder  in  integral  form  T{x\  p)  = 
A^t,  p)  +  fp  T^p\x]  f)  df.  We  can  now  bound  the  first  term  on  the  right  hand  side  of  (3)  by 

<  r  f|’;‘  <4) 

Jr-  U3-1!-  L=o(H)  P- 

for  all  p  G  V.  For  the  second  term  in  (3),  we  obtain 

I \APAt*iP)  -Pm^p)  l|Loo(n)  <  ||  APjt(t*,p)  -APjtm(t*,p)  ||ioo(n)  +  \\Ap^m{t*  ,  p)  -PM{",p)  ||ioo(fJ)  (5) 

for  all  p  G  V.  For  the  first  term  in  (5)  we  note  that 


>  P)  -  A'kST* ,  P) 


<  sup 
re* 


P~ 1  j 

ST'  Z3* 

si 

3=0  J' 


\/p  G  V. 


From  the  definition  of  the  characteristic  functions  ,  we  obtain  JZjZ o  (x>  t*)  ^~j\  ^  ~  p m(x ;  p)  = 

J2rn= l  ZZPjZo^M  (fm;  T*)  (M~y  )J  -  tFM(tm;p )  V^(:c).  We  then  invoke  the  interpolation  property  (for 

any  non-negative  integer  j)  (tm-,  p)  =  P^\tm;  p),  1  <  m  <  M,  and  the  definition  of  the  Lebesgue 
constant  AM,  to  bound  the  second  term  in  (5)  by 

ll-^M(r*’^)  -  11^,00(0)  <  \\-ApAtZp) -•iF(-;/3)llLoo(n)  am  <  \/pev.  (7) 

The  desired  result  (for  P  =  1)  directly  follows. 


We  make  several  remarks  concerning  this  result.  First,  we  may  choose  p  such  that  the  two  terms  in 
(2)  balance — a  higher  p  will  reduce  the  contribution  of  the  first  term  but  will  increase  the  contribution 
of  the  second  term.  Second,  we  note  that  the  bound  8m, P  is  ^.-independent.  We  can  readily  develop  a 
p-dependent  bound  by  replacing  p$  with  the  actual  distance  between  p  and  the  closest  r  G  d>;  this 
/i-dependent  bound  can  serve  (i)  to  adaptively  construct  an  economical  point  set  d>,  and  (ii)  to  replace 
the  true  (expensive)  error  in  the  greedy  identification  of  the  EIM  spaces  W^.  Third,  we  can  increase 
the  sharpness  of  the  bound  by  localizing  the  derivative  bounds  ap:  this  is  best  achieved  through  an  “/ip” 
approach  for  the  EIM;  we  note  that  the  “/ip”  framework  developed  in  [2]  for  the  reduced  basis  method 
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Figure  1.  Error  bounds  5m, p  for  P  =  1  and  p  =  1,2,  3, 4  with  =  41  (left)  and  =  141  (right). 

readily  adapts  to  the  EIM  (see  also  [5]  for  an  alternative  approach).  Fourth,  we  note  that  in  the  “limit” 
/?$  — >  0  the  effectivity  of  the  bound  approaches  unity;  of  course,  we  will  never  in  practice  let  — >  0 
because  this  implies  the  computation  of  the  interpolant  at  every  point  in  T>.  Fifth,  we  note  that  our  bound 
at  no  point  requires  computation  of  spatial  derivatives  of  the  function  to  be  approximated. 

We  conclude  this  section  by  summarizing  the  computational  cost  associated  with  8m,p-  We  assume 
that  the  bounds  ap  can  be  obtained  analytically.  We  compute  A m  in  0{M2J\f )  operations,  and  we 
compute  the  interpolation  errors  || —  ^m\'>  T)l|z/00(n)>  0  <  \P\  <  p  —  1,  for  all  r  G  <F,  in 
0(n$MAf)  Ej=o  card(Adj3)  operations  (we  assume  M  <C  W);  certainly  the  growth  of  Mp  will  preclude 
large  P.  Note  the  computational  cost  is  “offline”  only — the  bound  is  valid  for  all 


3.  Numerical  Results 

We  shall  consider  the  empirical  interpolation  of  a  Gaussian  function  p)  over  two  different  parameter 
domains  D  =  V\  and  V  =  V n.  The  spatial  domain  is  O  =  [0,1]2;  we  introduce  a  triangulation 
with  A f  =  2601  vertices.  We  shall  compare  our  bound  with  the  true  interpolation  error  over  the  parameter 
domain.  To  this  end,  we  define  the  maximum  error  Em  =  max^gHtrain  eM{p)  and  the  average  effectivity 
f)M,P  =  mean^gEtest 8m,p/ &m (m) 5  here,  eM{p)  =  p)  ~  p)\\l°°(Q),  and  Htest  C  V  is  a  test  sample 

of  finite  size  nntest. 

We  first  consider  the  case  V  =  V\  =  [0.1, 1]  and  hence  P  =  1;  we  let  T{x\  /i)  =  T\[x\  fi)  =  exp(— ((®i  — 
0.5)2  +  (x2  —  0.5)2)/2 /x2).  We  introduce  an  equidistant  train  sample  Strain  C  V  of  size  500;  we  take 
Hi  =  1  and  pursue  the  EIM  with  Mmax  =  12.  In  Figure  1  we  report  Em  and  8m, P >  P  =  1,2, 3, 4,  for 
1  <  M  <  Afmax;  we  consider  =  41  and  =  141  (p$  =  1.125 E-2  and  w  3.21  E-3,  respectively). 
We  observe  that  the  error  bounds  initially  decrease,  but  then  “plateau”  in  M.  The  bounds  are  very  sharp 
for  sufficiently  small  M,  but  eventually  the  first  term  in  (2)  dominates  and  compromises  the  sharpness 
of  the  bounds;  for  larger  p,  the  bound  is  better  for  a  larger  range  of  M.  We  find  that  1  <  A m  <  5.18  for 
1  <  M  <  Mmax  and,  for  the  case  p  =  4  with  n«j>  =  141,  T)m,p  ~  0(10)  (n=test  =  150)  except  for  large  M. 
The  modest  growth  of  the  Lebesgue  constant  is  crucial  to  the  good  effectivity. 

We  next  consider  the  case  V  =  Vn  =  [0.4,  0.6]2  and  hence  P  =  2;  we  introduce  T  =  T\\[x\p)  = 
exp(— ((xi  —  P(i))2  +  (%2  —  M(2})2)/2(0.1)2),  where  p  =  (/t(i), P(2))-  We  introduce  a  deterministic  grid 
Strain  CPof  size  1600;  we  take  pi  =  (0.4,  0.4)  and  pursue  the  EIM  with  Mmax  =  60.  In  Figure  2  we  report 
Em  and  Sm,p ,  P  =  1,2, 3, 4,  for  1  <  M  <  Afmax;  we  consider  =  100  and  =  1600  (p$  «  1.57E-2  and 
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Figure  2.  Error  bounds  8m, p  for  P  =  2  and  p  =  1,  2,  3,4  with  n$  =  100  (left)  and  =  1600  (right). 

3.63E-3,  respectively).  We  observe  the  same  behavior  as  for  the  P  =  1  case:  the  errors  initially  decrease, 
but  then  “plateau”  in  M  depending  on  the  particular  value  of  p.  We  find  that  1  <  Am  <  39.9  and,  for 
the  case  p  =  4  with  =  1600,  f]M,p  ~  0(10)  (n ste8t  =  225)  for  1  <  M  <  Mmax. 

Our  results  demonstrate  that  we  can  gainfully  increase  p — the  number  of  terms  in  the  Taylor  series 
expansion — in  order  to  reduce  the  role  of  the  first  term  of  Sm,p  and  to  limit  the  size  of  '!>.  We  also  note 
that  for  the  examples  presented  here  the  terms  in  the  sum  of  (2)  are  well  behaved,  even  though  (for 
our  P  =  2  example  in  particular)  it  is  not  obvious  that  the  space  W m  contains  good  interpolants  of  the 
functions  \(3\  0. 
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